############## 
############## 
############## Models
remove(list = ls())

base::library(conflicted)
base::library(tidyverse)
conflict_prefer("filter","dplyr")
base::library(ggplot2)
base::library(dplyr)
base::library(here)
conflict_prefer("here", "here")
base::library(systemfit)
base::library(stargazer)
base::library(DataCombine)

load(here("Data", "Model_Data.RData"))
glimpse(final_data)

summary(final_data$death_rate)
table(final_data$year)
hist(final_data$death_rate)

final_data <- final_data %>% 
  mutate(time_trend = ifelse(year == 2004, 1, NA),
         time_trend = ifelse(year == 2008, 2, time_trend),
         time_trend = ifelse(year == 2012, 3, time_trend),
         time_trend = ifelse(year == 2016, 4, time_trend),
         time_trend = ifelse(year == 2020, 5, time_trend)) 

library(fixest)
#?feols

ols_model1 <- lm(rep_share2 ~ death_rate + unemp_rate + log(defl_pcincome) +
                   share_black + share_hisp + share_asian + share_young + share_elder + 
                   log1p(share_college) + log1p(share_manuf) +
                   log(tot_pop) + as.factor(rural_urban) + state + state_partiesdiff + 
                   rep_inc + incumbent_cand + time_trend + lag_rep_share2,
                final_data); summary(ols_model1)


fe_model <- feols(rep_share2 ~ death_rate + unemp_rate + log(defl_pcincome) +
                    share_black + share_hisp + share_asian + share_young + share_elder + 
                    log1p(share_college) + log1p(share_manuf) +
                    log(tot_pop) + as.factor(rural_urban) + state + state_partiesdiff + 
                    rep_inc + incumbent_cand + time_trend + lag_rep_share2 | county,
                  data = final_data); summary(fe_model)


base::library(stargazer)

summary(ols_model1)

library(sandwich)
ols_model2 <- coeftest(ols_model1, vcov = vcovCL, cluster = ~ county)
ols_model2


stargazer(ols_model1, ols_model2, ols_model1,
          type = "latex",
          title = "OLS and Fixed Effect Estimates",
          dep.var.labels = c("Republican Vote Share", "Republican Vote Share", "Republican Vote Share"), 
          column.labels = c("State Fixed Effecs", 
                            "Clustered SE by County", 
                            "County Fixed Effects"),
          no.space = FALSE, single.row = T,
          covariate.labels = c("Overdose Death Rate",
                               "Local Unemployment Rate",
                               "Log PC Income",
                               "Share Black",
                               "Share Hispanic",
                               "Share Asian",
                               "Share Young Voters",
                               "Share Elder Voters",
                               "Log Share College Degree",
                               "Log Share Manufacture",
                               "Log Population",
                               "Close Election (State)",
                               "Republican Presidency",
                               "Incumbent Candidate",
                               "Lagged Republican Vote Share",
                               "Constant"),
          omit = c("rural_urban", "state"),
          omit.labels = c("Rural-Urban Code",
                          "State Fixed Effects"))



